# 导入必要的库
import numpy as np  # 用于数值计算和数组操作
import matplotlib.pyplot as plt  # 用于绘图和可视化
from sklearn.datasets import load_iris  # 加载Iris数据集
from sklearn.preprocessing import StandardScaler  # 数据标准化工具
from sklearn.metrics import silhouette_score, adjusted_rand_score  # 聚类评估指标
from scipy.cluster.hierarchy import dendrogram, linkage, fcluster  # 层次聚类工具
from sklearn.model_selection import train_test_split  # 划分训练集和测试集
from sklearn.decomposition import PCA  # 主成分分析，用于降维

# 设置matplotlib支持中文（解决中文显示问题）
plt.rcParams['font.sans-serif'] = ['SimHei']  # 使用黑体显示中文
plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示问题，避免乱码

# 1. 数据准备
# 加载Iris数据集，包含150个样本，4个特征（萼片长度、宽度，花瓣长度、宽度）
iris = load_iris()
X = iris.data  # 特征矩阵：150 x 4，包含所有样本的特征数据
y = iris.target  # 真实标签：0, 1, 2（对应setosa、versicolor、virginica），仅用于评估

# 数据清洗：检查并处理缺失值
if np.any(np.isnan(X)):  # 检查数据中是否包含NaN值
    print("警告：数据包含缺失值，进行均值填充...")
    # 按列计算均值，避免空数组问题
    column_means = np.nanmean(X, axis=0)
    # 如果某列全为NaN，用0填充
    column_means = np.where(np.isnan(column_means), 0, column_means)
    # 用列均值填充NaN值，保持数据形状
    X = np.where(np.isnan(X), np.tile(column_means, (X.shape[0], 1)), X)
else:
    print("数据无缺失值，跳过清洗步骤")

# 数据标准化：使用Z-score标准化，消除不同特征的量纲影响
scaler = StandardScaler()  # 初始化标准化器，计算均值和标准差
X_scaled = scaler.fit_transform(X)  # 标准化数据，均值为0，方差为1

# 划分训练集和测试集（仅用于评估聚类效果）
X_train, X_test, y_train, y_test = train_test_split(X_scaled, y, test_size=0.2, random_state=42)

# 2. 计算距离矩阵
# 使用scipy的linkage函数，自动计算距离矩阵并进行初始层次聚类
# 尝试多种链接方法以选择最佳聚类效果
methods = ['single', 'complete', 'average', 'ward']  # 不同的链接方法
best_ari = -1  # 初始化最佳调整Rand指数
best_method = None  # 初始化最佳方法
best_y_pred = None  # 初始化最佳聚类标签
best_Z = None  # 初始化最佳链接矩阵

for method in methods:
    # linkage函数自动计算距离矩阵并进行层次聚类
    # method参数指定链接方法：single(单链接)、complete(全链接)、average(平均链接)、ward(方差最小化)
    Z = linkage(X_scaled, method=method)
    # 3. 初始化聚类：linkage隐式初始化每个样本为一个簇，无需额外步骤

    # 4. 合并聚类与5. 重复合并过程
    # linkage函数逐步合并最近的簇，Z记录所有合并信息
    # 选择聚类数：根据Iris数据集的真实类别数，设为3
    n_clusters = 3
    y_pred = fcluster(Z, t=n_clusters, criterion='maxclust')  # 切割树形图生成聚类标签

    # 评估聚类效果，使用调整后的Rand指数
    ari = adjusted_rand_score(y, y_pred)  # 比较聚类标签与真实标签的相似度
    print(f"链接方法: {method}, 调整后的Rand指数: {ari:.4f}")

    # 记录最佳方法和结果
    if ari > best_ari:
        best_ari = ari
        best_method = method
        best_y_pred = y_pred
        best_Z = Z

print(f"\n最佳链接方法: {best_method}, 最佳调整Rand指数: {best_ari:.4f}")

# 6. 生成树形图（基于最佳方法）
# 绘制层次聚类的树形图，展示合并过程和距离变化
plt.figure(figsize=(12, 6))
dendrogram(best_Z, labels=[str(i) for i in range(len(X))], leaf_rotation=90)
plt.title(f'层次聚类树形图 - 链接方法: {best_method}')
plt.xlabel('样本索引')
plt.ylabel('距离')
plt.show()

# 7. 结果分析与评估
# 计算轮廓系数，评估簇内紧密度和簇间分离度
silhouette_avg = silhouette_score(X_scaled, best_y_pred)
print(f"轮廓系数 (Silhouette Score): {silhouette_avg:.4f}")

# 输出调整后的Rand指数（已通过最佳方法计算）
print(f"调整后的Rand指数 (Adjusted Rand Score): {best_ari:.4f}")

# 5. 可视化聚类结果
# 使用PCA将4维数据降维到2D，以便可视化
pca = PCA(n_components=2)  # 初始化PCA，降维到2维
X_pca = pca.fit_transform(X_scaled)  # 转换数据到2D空间

# 绘制聚类结果散点图
plt.figure(figsize=(10, 6))
scatter = plt.scatter(X_pca[:, 0], X_pca[:, 1], c=best_y_pred, cmap='viridis', edgecolor='k')
plt.title(f'层次聚类结果 (簇数={n_clusters}, 链接方法: {best_method})')
plt.xlabel('第一主成分')
plt.ylabel('第二主成分')
plt.colorbar(scatter, label='簇标签')
plt.show()

# 绘制真实标签对比散点图
plt.figure(figsize=(10, 6))
scatter = plt.scatter(X_pca[:, 0], X_pca[:, 1], c=y, cmap='viridis', edgecolor='k')
plt.title('真实标签分布')
plt.xlabel('第一主成分')
plt.ylabel('第二主成分')
plt.colorbar(scatter, label='真实类别')
plt.show()

# 8. 应用与扩展
# 统计每个簇的样本数量
cluster_counts = np.bincount(best_y_pred)
print("\n每个簇的样本数量:")
for cluster in range(max(best_y_pred)):  # 仅遍历存在的簇
    count = cluster_counts[cluster]
    if count > 0:  # 仅打印非空簇
        print(f"簇 {cluster} 包含 {count} 个样本")

# 计算每个簇的平均特征，作为后续分析的基础
# 过滤空簇，确保计算安全
cluster_centers = np.array([np.mean(X_scaled[best_y_pred == i], axis=0)
                           for i in range(n_clusters) if np.any(best_y_pred == i)])
print("\n每个簇的平均特征（标准化后）:")
for i, center in enumerate(cluster_centers):
    print(f"簇 {i} 平均特征: {center}")

# 扩展应用示例：输出簇的特征差异
print("\n簇特征差异分析（簇0与簇1的差值）：")
if len(cluster_centers) >= 2:
    diff = cluster_centers[0] - cluster_centers[1]
    print(f"簇0与簇1特征差异: {diff}")
else:
    print("簇数不足，无法计算差异")